Introduction

Here we will recreate the ZooMSS model (version 2) in Heneghan et al. (2020 in preparation) using mizer.

We begin with some setup of required packages.

#get required packages
library(devtools)
## Warning: package 'devtools' was built under R version 3.6.3
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 3.6.3
#most up to date master branch of mizer
#install_github("sizespectrum/mizer")
#install_github("astaaudzi/mizer-rewiring", ref = "temp-model-comp")
#documentation here:
#https://sizespectrum.org/mizer/dev/index.html
library(mizer)
require(tidyverse)
## Loading required package: tidyverse
## Warning: package 'tidyverse' was built under R version 3.6.3
## -- Attaching packages --------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.0
## v tidyr   1.1.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 3.6.3
## Warning: package 'tibble' was built under R version 3.6.3
## Warning: package 'tidyr' was built under R version 3.6.3
## Warning: package 'readr' was built under R version 3.6.3
## Warning: package 'purrr' was built under R version 3.6.3
## Warning: package 'dplyr' was built under R version 3.6.3
## Warning: package 'stringr' was built under R version 3.6.3
## Warning: package 'forcats' was built under R version 3.6.3
## -- Conflicts ------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
remotes::install_github("sizespectrum/mizerExperimental")
## Skipping install of 'mizerExperimental' from a github remote, the SHA1 (2e116191) has not changed since last install.
##   Use `force = TRUE` to force installation
library(mizerExperimental)

Data

Next, let’s load in and check the data used in the ZooMSS model. TODO: add in some biomass data

#read in group data
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
#groups<-merge(groups,dat,by.x="species",by.y="group",all=T)

# have a look at plot
# plot(groups$w_inf,groups$biomass.tperkm2,xlab="log Maximum Weight [g]", ylab=" log Total Biomass", log="xy",col="blue",pch=16)
# text(groups$w_inf,groups$biomass.tperkm2,labels=groups$species,cex=0.5)

# could plot the paramter allometries here to explore

Set-up mizer model

Next let’s read in the parameters from ZooMSS.

#read groups again for southern ocean model, this time subsetting key groups
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min #convert from log10 values
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
groups$h <- 100000 #maximum feeding level - very large value to reflect the unlimited feeding implicitly allowed in ZooMSS

#groups <- readRDS("data/groups.RDS")[-1,]
#check fails these tests:
groups$w_mat25 >= groups$w_mat #note we don't have w_mat25. Not sure if it's necessary/desirable. Probably not.
## logical(0)
groups$w_mat25 <= groups$w_min
## logical(0)
groups$w_inf <= groups$w_mat25
## logical(0)
#[-1,]
#fix one value 
#groups[7,"w_min"]<-4000

# read interaction matrix
# get the interaction matrix - actually I think we can leave this out. Default is all 1s, which is the same as in ZooMSS. Included for completeness; it may be useful in future to keep this in.
theta <- readRDS("data/zoomss_inter.rds")[,-1]
#[-1,-1]

We will pass these parameters to mizer to set up a new multispecies model.

TODO: adjust parameters here.

mf.params <- newMultispeciesParams(species_params=groups,
                                   interaction=NULL, #NULL sets all to 1
                                   no_w = 178, #number of zoo+fish size classes;
                                   kappa = 1e15, #set huge to recreate unbounded feeding on phyto in ZooMSS
                                   min_w_pp = 10^(-14.5), #minimum phyto size
                                   w_pp_cutoff = 10^(-8), #maximum phyto size
                                   r_pp = 20, #Coefficient of the intrinsic resource birth rate, set large for now
                                   n = 0.7, #The allometric growth exponent used in ZooMSS
                                   z0pre = 1, #external mortality
                                   z0exp = 0.3,
                                   resource_dynamics = "resource_constant",
                                   #RDD = constantRDD(species_params = groups) #first go at this
                                   #pred_kernel = ... #probably easiest to just import this/pre-calculate it, once dimensions are worked out
)
## Note: No ks column so calculating from critical feeding level.
## Note: Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
# # mf.params@species_params$metab[1] <-0.01
# mf.params@species_params$h[] <- 10*mf.params@species_params$ks[] 
# 
# # mf.params@species_params$gamma[1] <- 1000*(mf.params@species_params$gamma[1])
# mf.params@species_params$alpha[] <- 0.6
# #mf.params@species_params$z0[1] <- 0.01
# 
# #higher mean betas
#  mf.params@species_params$beta[mf.params@species_params$species == "large.divers"] <-6000
#  mf.params@species_params$beta[mf.params@species_params$species == "apex.predators"] <-1000
# # # wider feeding kernels
#  mf.params@species_params$sigma[] <- 2
#  mf.params@species_params$sigma[mf.params@species_params$species == "baleen.whales"] <- 4
# 
# 
# # mf.params@species_params$w_mat25[1] <- 0.005*mf.params@species_params$w_inf[1]
# mf.params@species_params$erepro <- ifelse(mf.params@species_params$erepro==1,0.001,mf.params@species_params$erepro)
# 
# # mf.params@interaction[] <-0.1
# #mf.params@interaction[,1] <-0
# #mf.params@interaction[1,1] <-1
# 
# #mf.params@species_params$interaction_p[2:11] <-0.5
# 
# # hand tuning to gte in line with data
# mf.params@species_params$R_max <- 10*mf.params@species_params$R_max
# 
# 
# # needs to be very different for marine mammals? use different repro assumptions (fixed offspring density per year)?
# 
# 
# 
# # mf.params@species_params$w_mat <- 0.5*mf.params@species_params$w_inf
# # mf.params@species_params$w_mat[mf.params@species_params$species == "flying.birds"] <- 0.75*mf.params@species_params$w_inf[mf.params@species_params$species == "flying.birds"]
# # mf.params@species_params$w_mat[mf.params@species_params$species == "small.divers"] <- 0.75*mf.params@species_params$w_inf[mf.params@species_params$species == "small.divers"]
# # mf.params@species_params$w_mat[mf.params@species_params$species == "medium.divers"] <- 0.75*mf.params@species_params$w_inf[mf.params@species_params$species == "medium.divers"]
# # mf.params@species_params$w_mat[10:12] <- 0.75*mf.params@species_params$w_inf[10:12]
# # setParams(mf.params)
# 
# mf.params<- setParams(mf.params,kappa=1e5) # take the new paramewters and change kappa


sim <- project(mf.params, t_max=500,dt = 0.1)
plot(sim)

#plotlyGrowthCurves(sim,species="macrozooplankton")
plotlyFeedingLevel(sim)
# feeding level satioation for some groups, except for the seabirds
# macrozooplankton - they are not growing enough,why?
#tuneParams(mf.params)
plotlyGrowthCurves(sim,percentage = T)
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 0.127549
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 0.00198096
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps

## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 0.00198096
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps

## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
plotlySpectra(sim)

—–Old stuff only below this line—–